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Abstract 



We present Monte Carlo simulation results on the equilibrium relaxation dy- 
namics in the two dimensional lattice Coulomb gas, where finite fraction / of 
the lattice sites are occupied by positive charges. In the case of high order 
rational values of / close to the irrational number 1 — g (g = (\/5 — l)/2 is 
the golden mean), we find that the system exhibits, for wide range of tem- 
peratures above the first-order transition, a glassy behavior resembling the 
primary relaxation of supercooled liquids. Single particle diffusion and struc- 
tural relaxation show that there exists a breakdown of proportionality between 
the time scale of diffusion and that of structural relaxation analogous to the 
violation of the Stokes-Einstein relation in supercooled liquids. Suitably de- 
fined dynamic cooperativity is calculated to exhibit the characteristic nature 
of dynamic heterogeneity present in the system. 

PACS No.: 64.70.Pf, 64.60.Cn 
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I. INTRODUCTION 

The dynamics of supercooled liquids ap- 
proaching the glass transition remains one of 
the most fundamental problems in condensed 
matter physics 0. Some of the prominent 
dynamic features in supercooled liquids in- 
clude the enormous increase in relaxation 
time scale with lowering temperature, and 
the nonexponential relaxation in the response 
to an external perturbation. In addition 
to these features, an anomaly in transport 
properties such as breakdown of the Stokes- 



Einstein (SE) relation in highly supercooled 
liquids has been observed in experiments |2j 
and simulations @^5|] . Although there exist 
some theoretical attempts [p|-HT0|], the under- 
lying microscopic mechanism for the viola- 
tion of the SE relation is not well understood. 
Recently, there have been many experimen- 
tal and simulational studies of supercooled 
liquids that demonstrate the existence of ki- 
netic heterogeneity which was often invoked 
to explain the origin of the non-exponential 
relaxation as well as the breakdown of the SE 



relation 11 . 



In relation to these questions on micro- 
scopic slow dynamic features in supercooled 
liquids, we deemed it worthwhile to investi- 
gate whether similar dynamic features can be 
found in simpler lattice spin systems or lat- 
tice gas systems. In this work, we show that 
the aforementioned features of supercooled 
liquids, i.e., slowing-down, non-exponential 
relaxation and the (analogue) of the break- 
down of the SE relation, are also observed 
in a two-dimensional (2D) lattice Coulomb 
gas (LCG) system. We also find that the re- 
laxation of the system exhibits a spontaneous 
appearance of spatial heterogeneity, which we 
argue is the underlying cause for the non- 
exponential relaxation and the breakdown of 
the SE relation. 

In recent years, there have been some ef- 
forts to find glassy dynamic features in the 
lattice spin systems with nonrandom inter- 
actions 111211 . One of the well-known exam- 



ples of disorder-free lattice model system is 
uniformly frustrated XY (UFXY) models in 
two dimensions, which serve as a model for 
two dimensional arrays of Josephson junc- 
tions under the influence of uniform trans- 
verse magnetic fields. Recent work |H| has 
shown that, irrespective of the true nature 
of the low temperature phase of this system, 
the equilibrium dynamics of UFXY model in 
the intermediate range of the temperature 
for frustration parameter / near 1 — g = 
(3 — Vo)/2 — 0.382, exhibits a close anal- 
ogy to that of supercooled liquids. Both spin 
and vortex dynamics show stretched expo- 
nential relaxations with temperature depen- 
dent stretched exponents. In order to inves- 
tigate the dynamics of this system in more 
detail, we attempted to calculate the self dif- 
fusion properties of vortices. However, it 
turned out to be numerically ambiguous and 
tricky to trace the trajectories of individual 
vortices. This is beacause individual vor- 
tex around a plaquette is defined in terms of 
phase angles and one probes the movement 



of individual vortices not directly but only 
indirectly through changes of phases, which, 
at times especially when multi- vortex motion 
occurs, makes it ambiguous to determine the 
original position of a vortex corresponding to 
a new neighboring vortex. 

One way to overcome this difficulty was to 
map the UFXY model onto a LCG via Vil- 
li! 



lain transformation 



where the positive 



charges in the LCG correspond to the posi- 
tive current vortices in UFXY models. One 
can readily probe the diffusive dynamics of 
charges without ambiguity in the LCG un- 
like the case of UFXY model. Hence we can 
investigate both the structural relaxation dy- 
namics and self diffusion dynamics of individ- 
ual vortices in LCG with relative ease. 

With this advantage, we have numerically 
investigated the equilibrium relaxation dy- 
namics and diffusion characteristics of LCG 
with charge density factor / near 1 — g ~ 
0.382. We observe that for some range of 
temperatures above the first order transition, 
the equilibrium relaxation exhibits slow dy- 
namic features such as stretched exponential 
relaxation and breakdown of proportionality 
between diffusive time scale and structural 
relaxation time scale. 

It was a common belief that the 2D 
UFXY model and the corresponding LCG be- 
long to the same universality class with essen- 
tially the same phase transition properties, 
ground state symmetry, for example. How- 
ever, recent work on LCG by Gupta, Tei- 
tel, and Gingras (GTG) |L5| and also an- 
other work on UFXY model by Denniston 
and Tang (DT) [Jl6| , showed that there ex- 
ist some difference between the two model 
systems especially in the case of dense frus- 
tration. Both model systems exhibit first or- 
der transition but the low temperature vortex 
configurations in UFXY models are different 
from the charge configurations of the corre- 
sponding LCG for / near 1 — g ~ 0.382. The 
underlying cause for this breakdown of Vil- 
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lain approximation in the limit of dense frus- 
tration is not known, but probably it is re- 
lated to application of spin-wave integration 
to systems having many metastable states 
with similar energies, that may cause neglect 
of multi- vortex correlations. 

Special interest has been given to the case 
of / approaching 1 — g |T7|JI8| . Consider 
a system where / equals po/<?o (Po an d <?o 
are relative primes) which is a rational ap- 
proximant to 1 — g. Here, in the case of 
UFXY model, DT argues that the low tem- 
perature vortex configuration has lattice pe- 
riodicity which is of order q^, i. e., much larger 
than q . On the other hand, in the case 



of LCG, GTG [0 showed, via Monte Carlo 
(MC) simulations, that the low temperature 
charge configurations are characterized by ar- 
rangements of diagonal stripes that are either 
completely filled, completely empty, or par- 
tially filled with charges that are quite differ- 
ent from those vortex configurations in the 
corresponding UFXY model. However, GTG 
did not enumerate the exact patterns of low 
temperature charge configurations (such as 
spatial periodicity) for general cases of dense 
charge filling. In this work, we find that, for 
the values of / between 1/3 and 2/5, there 
exist a simple regularity in the low temper- 
ature charge configuration which consists of 
periodic arrangements of combinations of two 
out of three types of striped charge patterns 
(see Section III). 

For wide range of quenching tempera- 
tures above the first order transition T c , 
the equilibrium relaxation continues to slow 
down with lowering temperature, and the 
form of the relxations are characterized by 
the stretched exponential with temperature- 
dependent exponents. Moreover, we ob- 
serve that the model exhibits a separation 
of the two characteristic time scales, i.e., 
the time scale of single particle diffusion and 
that of structural relaxation. This feature 
is quite analogous to the breakdown of the 



SE relation observed in supercooled liquids. 
Stretched exponential relaxation is observed 
to be accompanied by interesting dynamic 
heterogeneity in the system. It appears that 
the kinetic heterogeneity is the underlying 
reason for both the stretched exponetial re- 
laxation and the separation of the relaxation 
and diffusion time scales. 

A convenient measure for dynamic het- 
erogeneity is the so called dynamic coopera- 
tivity [^] of the particle motions. This mea- 
sures reduction of the effective degrees of free- 
dom. One interesting result from our simula- 
tions is that the magnitude of the velocity (or 
displacement vector) exhibits strong increase 
in cooperativity of the particle motions. On 
the other hand, the displacement vector it- 
self shows cooperativity a little smaller than 
unity due to anti-correlations in the direc- 
tion of particle motions. This means that 
the system can be divided into highly mobile 
regions and relatively inert regions resulting 
in highly inhomogeneous local mobility dis- 
tribution. However, there is no macroscopic 
flow of particles that will generate long range 
positive correlations between the directions 
of flows of particles. 

When quenched to a temperature below 
T c , the system is always found to undergo 
phase ordering via slow coarsening processes. 
The system therefore does not remain in a 
supercooled state. Rather it becomes slowly 
crystalized. It should be emphasized that in 
this system it is the relaxation for the temper- 
atures above T c that exhibits slow dynamic 
behavior which shares some common features 
with that of supercooled liquids. 

II. MODEL AND SIMULATION 
METHODS 

General 2D LCG [^(J is described by the 



following Hamiltonian that can be mapped 
from UFXY model by means of Villain trans- 
formation 111411 , 
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ij 

where is the distance between the sites i 
and j, and the magnitude of charge Qi at 
site i can take either 1 — / or — / , where 
/ corresponds to the frustration parameter 
in the related XY models. Charge neutrality 
condition J2i Qi — implies that the number 
density of the positive charges is equal to /. 
As was mentioned above, we can thus view 
the system as a lattice gas of N ■ / charges of 
unit magnitude upon uniform negative back- 
ground charges of charge density — / (N = L 2 
is the total size of the system with the linear 
dimension L). The lattice Green's function 
G(rij) solves the equation 
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where A 2 is the discrete lattice Laplacian and 
A is the screening length which, in normal 
case of no screening, is set to an infinity. For 
the case of usual Villain transformation of 
UFXY model, we have A = oo. But it is in- 
cluded in this equation for generality. Since, 
in this work, we restrict our attention to only 
square lattice with periodic boundary condi- 
tions, G(r) is given by 
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where k are the allowed wave vectors with 
kfj, = (27rn M /L), with — 0, 1, . . . , L — 1. In 
the case of infinite screening length, for large 
separation r, one gets G(r) ~ — lnr pl |. In 
this work, we consider the limiting case of 
A — > oo only. 

In our MC simulations, the initial disor- 
dered random configuration is updated ac- 
cording to the standard Metropolis algorithm 
by selecting a positive charge at random and 
moving it over to one of the nearest neighbor 



(NN) or next nearest neighbor (NNN) sites 
fl5fl . We find that this NNN hopping al- 
gorithm is particularly effective in terms of 
simulation time as compared with NN hop- 
ping alone, as was emphasized in [|Tj| . More- 
over, at low temperature, NN hopping alone 
presented severe energy barriers to the mo- 
tions of charges in the case of relatively dense 
Coulomb gas, i.e., f approximately larger 
than 1/3. 

The presented results are averages over 
100 ~ 500 different random initial configura- 
tions depending on the temperature. In or- 
der to ensure that equilibration is achieved, 
we calculate the two-time charge density au- 
tocorrelation function and locate the waiting 
time beyond which the autocorrelation func- 
tion no longer depends on the waiting time. 
As for the values of the charge density pa- 
rameter /, we chose / = 55/144 ~ 0.3819 
that is close to / = 1 — g, and square lattices 
of linear size L = 36 is chosen with periodic 
boundary conditions. This value of / is cho- 
sen as a simple rational value that satisfies 
the two conditions of both being close to 1— g 
and being commensurate with the lattice pe- 
riodicity 12 as explained in Section III. We 
found that qualitative features of relaxation 
dynamics are the same for other nearby val- 
ues of the frustration /. 

III. SIMULATION RESULTS AND 
DISCUSSIONS 

A. First order transition and low 
temperature configuration 

We first discuss the equilibrium phase 
transition and charge configuration of the 
system. As was first shown by GTG, we 
also find that there exist a first order tran- 
sition in LCG with / near 1 — g. Figure 1 
shows temporal snapshots of charge config- 
urations evolving from disordered state into 
ordered configuration after being quenched to 
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a temperature T = 0.026. First order na- 
ture of the phase transition can easily be con- 
firmed by enumerating the histogram of en- 
ergies P(E) near the transition temperature 
22]]. P{E) is obtained by counting the oc- 



currences of energies for each of the equally 
spaced energy bins while performing the equi- 
librium Monte Carlo simulations (via simple 
Metropolis algorithm). For a system with 
first order transition, the energy histogram 
P(E) becomes bimodal near the transition 
temperature corresponding to a mixture of 
ordered state (with lower energy) and a dis- 
ordered state (with higher energy). The tran- 
sition temperature T c can be determined by 
locating the temperature where the subar- 
eas under the two peaks are equal. Fig- 
ure 2 shows two histograms near the tran- 
sition temperature, where we could estimate 
the transition temperature approximately as 
T c ~ 0.0316. Since we did not attempt a de- 
taled analysis (including a finite size scaling) 
of the histogram, we think that this estimate 
value of the transition temperature should 
not be considered too seriously for its pre- 
cision. 

We find empirically that there exist a sim- 
ple regularity in the low temperature charge 
configuration in LCG (Fig. 3). For the case of 
values of / in the range 1/3 < / < 2/5, it is 
found that the low temperature configuration 
becomes quasi-one-dimensional with periodic 
striped patterns. In the cases of / = 1/3 and 
/ = 2/5 the ground state configurations are 
identical to the low temperature vortex con- 
figurations in the UFXY model. However, 
for values of / in between 1/3 and 2/5, the 
low temperature patterns are found to be, un- 
like the case of corresponding UFXY model, 
consisting of periodic arrangements of com- 
binations of two out of three types of striped 
charge patterns as follows. 

First component pattern (type I pattern) 
is a sequence of three diagonals which are 
empty, filled, and empty respectively (that 



may be denoted by (010) in our notation 
where 1 refers to a filled diagonal and refers 
to an empty diagonal). In other words, it is 
a pattern with single isolated diagonal filled 
with charges, that is neighbored by empty 
diagonals on both sides. Repetition of this 
pattern alone produces the ground state con- 
figuration for the case of / = 1/3 with spatial 
periodicity three. 

Second component pattern (type II pat- 
tern) consists of a sequence of five diago- 
nals that are empty, filled, empty, filled, and 
empty respectively, or (01010) in our nota- 
tion. This may be termed as a double filled 
diagonal because two filled diagonals are po- 
sitioned in parallel at second neighbor. This 
forms the basis of the ground state configu- 
ration for the case of / = 2/5 with lattice 
periodicity five. 

Lastly, the third component pattern (type 
III pattern) consists of a sequence of seven 
diagonals that are sequentially empty, filled, 
empty, partially filled, empty, filled, and 
empty i.e., (OlOpOlO) in our notation where 
p refers to a partially filled diagonal where 
only part of the diagonal sites are occupied by 
positive charges. This pattern is essentially 
a partially filled diagonal enveloped by two 
filled diagonals on both sides at second neigh- 
bor diagonal position, which may be termed 
as a channel structure. This can form a basis 
with spatial lattice periodicity seven. 

We leave the detailed description of the 
low temperature charge patterns for the full 
range of / values between 1/3 and 2/5 to 
the forthcoming publication [|23|| . And we de- 
scribe in this work the low temperature or- 
dered patterns for values of / around 1 — g 
only. Near the value of the filling ratio / = 
1 — g ~ 0.382, we find that, among the three 
patterns above, only two types (type II and 
type III patterns) participate in the stable 
charge configurations with the resulting spa- 
tial lattice periodicity depending on the com- 
bination of the two component patterns. 
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We find that there exist a value / = f c ~ 
0.381 which separates two regimes with dis- 
tinct low temperature striped patterns. For 
values of / in the range 0.36 ^ / ^ 0.381, 
the stable striped patterns turn out to have 
periodicity l p — 7 which consists of simple 
repetitions of channel structures (type III 
pattern). Note that this periodicity seven 
refers to the periodicity of the filled diago- 
nals only (neglecting the true periodicity in- 
cluding the charge configurations within the 
partially filled diagonals). 

On the other hand, for values of / in the 
range 0.381 ^ / ^ 0.39, the stable configura- 
tion exhibits a periodicity l p = 12, which con- 
sists of double filled diagonals (type II) and 
channels (type III) alternatingly placed. As 
the value of / continuously increases within 
the two regimes (in the above), the system 
in the low temperature stable configuration 
simply adjusts itself by accomodating the ex- 
tra charges into the partially filled diagonal 
channels and thereby changing the charge fill- 
ing within the channels. The dividing value 
of / = f c — 0.381 between the two regimes 
appears to correspond to the value 8/21 in 
which case the partially filled diagonals have 
filling density exactly equal to 2/3. Our sim- 
ulations show that the filling density 2/3 in- 
side the partially filled diagonal plays as a 
stability limit for the channel structures. Be- 
yond this limit, electrostatic instability prob- 
ably begins to set in, and rearrangement of 
the whole charge configuration occurs in or- 
der to form a new stable ordered patterns. 
As was also argued by GTG, in general, at 
much lower temperature T p (below T c ) the 
charges within the partially filled channels 
are expected to exhibit ordering, which would 
depend sensitively on rationality of the ex- 
act filling ratio of charges inside the partially 
filled diagonals. 

An important aspect of our simulations 
is that one has to choose the lattice size ap- 
propriately in order to match the periodicity 



of the true low temperature configuration in 
the thermodynamic limit. If, otherwise, one 
chooses a lattice size that is incommensurate 
with the periodicity (of striped patterns), 
then one ends up with defective charge con- 
figurations with patches of local ground state 
configurations. We think that this is prob- 
ably why GTG got two different equilibrium 
configurations when two different lattice sizes 
L = 26 and L = 52 are used for / = 5/13 
since these L's turn out to be incommensu- 
rate with the true periodicity l p = 12. 

When the screening length A is finite, then 
we find the low temperature configuration be- 
comes different from the case of no screening 
(A — > oo) in such a way that the partially 
filled diagonals gets rarer. The influence of 
the screening effect on the statics and the re- 
laxation dynamics needs further study. 

B. Equilibrium relaxation dynamics 

We now discuss the equilibrium relaxation 
dynamics of the model above first order tran- 
sition. In order to probe the structural re- 
laxation of charges, we measured the on-site 
charge autocorrelation functions, 

C(t) = (J2Q l (0)Q i (t))/N, (4) 
i=i 

where the bracket < • • • > represents an av- 
erage over different random initial configura- 
tions. 

Shown in Fig. 4a is the on-site charge au- 
tocorrelation function C(t) for temperatures 
from T = 0.1 down to T = 0.033. From 
this figure, we observe a slowing down in 
the structural relaxation for this temperature 
range. One can extract a characteristic time 
scale t(T) which, for example, is defined as 
C(t — t(T)) = 1/e for each temperature T. 
As Fig. 4b clearly shows, the temperature 
dependence of the relaxation time exhibits 
a non-Arrhenius behavior. We also checked 
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whether the so-called time-temperature su- 
perposition holds for the above autocorrela- 
tion functions, which is shown in Fig. 4c. We 
clearly see that time-temperature superposi- 
tion is systematically broken by the autocor- 
relation functions. This is consistent with the 
fact that the stretched exponents have depen- 
dence on temperature as is shown just below. 

We find that the relaxation pattern of 
the correlation function C(t) can be char- 
acterized by a power law relaxation C(t) = 
1 — At b ^ (known as the von Schweidler 
relaxation) in the early time regime and 
a stretched exponential relaxation C(t) = 
Cq(T) exp(— A't^^) in the late time regime. 
However, as the temperature gets higher, the 
regime of validity for early time power law 
relaxation was significantly reduced and we 
could better fit the early time relaxation with 
another stretched exponential form C(t) = 
exp(—A"t b ^). Of course for low tempera- 
ture regime, we could get b(T) ~ b'(T). 

Fig. 4d shows the temperature depen- 
dence of the fitted exponents. We see that 
non-exponentiality increases as the temper- 
ature decreases. These results clearly indi- 
cate that the equilibrium relaxation in the 2D 
LCG above T c closely resembles the primary 
relaxation of typical fragile liquids. 

One of the main characteristic features of 
the single particle dynamics is described by 
the mean square displacement ((Ar) 2 ), which 
is defined as 

((Arf) = (j:m)~r J (0)) 2 )/N Q , (5) 

where fj(t) is the position vector of the j- 
th charge at time t and Nq the total num- 
ber of charges. Figure 5 shows ((Ar) 2 ) for 
various temperatures. It exhibits an early 
time subdiffusive regime and crosses over into 
late time diffusive regime. Early time sub- 
diffusive behavior is thought to be coming 
from local frustrated motions of charges be- 
fore reaching an average displacement of unit 



lattice spacing. To test the proportionality 
of the two time scales, the structural relax- 
ation time scale r and the diffusion time scale 
D~ x , we plot the temperature dependence of 
the product ADt in Fig. 6. Here, we clearly 
see that the breakdown of the proportionality 
between the two time scales is observed for 
wide range of temperatures below T = 0.1 
and becomes stronger as the temperature is 
lowered. This separation of the two time 
scales is due to the weaker temperature de- 
pendence of the diffusion coefficient. That is, 
diffusion is relatively enhanced at lower tem- 
peratures. This is quite analogous to the vi- 
olation of the SE relation (D = T/arj, where 
a is a molecular length and t] is the viscos- 
ity of the liquid) observed in experiments on 
supercooled liquids 0. Here, we mention 
that there exists a correlation between the 
increase of non-exponentiality (as the tem- 
perature is lowered) and the increase of the 



product ADt at low temperatures [24 



If we suppose that there exists a sin- 
gle dominant relaxation mode in the sys- 
tem (and hence one relaxation time scale r), 
then we would obtain a simple exponential 
behavior for the relaxation function C(t) ~ 
e~ l l T . On the other hand, if the system con- 
sists of many regions with different relaxation 
times, then the relaxation function would be 
roughly some superposition of exponentials 
with a broad distribution of relaxation times, 
which would be in general not expressible in 
a simple exponential form, but in stretched 
exponetial form or other more complicated 
forms. 

The fact that there exists a breakdown 
of proportionality between r and D^ 1 can 
be interpreted in the following way that in- 
vokes dynamic heterogeneity. As the temper- 
ature is lowered, the system consists of many 
regions with different relaxation time that 
comes from different local mobilities. We 
can easily see that the structural relaxation 
time is dominated by the least mobile regions, 
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that is, by the regions with the longest re- 
laxation time. On the contrary, the average 
(long time) diffusion characteristics is domi- 
nated by the most mobile regions. In other 
words, the structural relaxation function and 
the self-diffusion function, respectively, are 
probing more or less opposite aspects of the 
relaxation behavior of the system. For an 
extreme example, one can imagine a system 
where half of the whole system is frozen (no 
motion of the component particles) while the 
remaining half of the system has finite relax- 
ation time with uniformly distributed mobile 
particles. For this system, the structural re- 
laxation time t would be infinite due to the 
frozen half of system, but inverse of the av- 
erage diffusion constant D~ x is finite due to 
the mobile part of the system, leading to an 
extreme breakdown of SE relation. Above 
simulation result, thus, can be interpreted as 
an evidence pointing toward the existence of 
a kinetic heterogeneity in the relaxation dy- 
namics and mobility of the system. 

In fact, the kinetic heterogeneity can be 
visualized in our system. Typical charge 
configuration at T = 0.033, as shown in 
Fig. 7, exhibits local striped patterns (or- 
dered domains) and interfacial regions due 
to mismatch between adjacent domains. For 
a fixed quenching temperature, the average 
size of these local domains reaches a certain 
length scale when the system equilibrates. 
After equilibration, the system structurally 
rearranges itself going from one configura- 
tion to another with local domains of sim- 
ilar length scale. Figure 8 shows the tra- 
jectories of moving positive charges over a 
time interval of 500 MC steps for T = 0.033 
(corresponding to Fig. 7). We can see that 
there exist local regions with actively moving 
charges and other regions with relatively im- 
mobile charges. Among the active regions, we 
can find those charges moving along partially 
filled diagonal channels. We also find some 
extended interfacial regions where no dis- 



cernible local order can be identified, that ex- 
hibit relatively high mobility. Enhancement 
of particle diffusion is probably due to the 
motions of charges along the partially filled 
diagonals as well as those fluidized motions in 
the extended interfacial regions. These fastly 
moving regions in surroundings of very slowly 
moving regions offer a specific example for 
spatial heterogeneity in glassy systems ||||, 
which was often thought of as the physical 
mechanism for breakdown of the SE relation. 

One simple way to quantify the degree of 
dynamic heterogeneity directly from the local 
motions of particles is to calculate the dy- 



namic cooperativity |19j for one particle dy- 
namic quantities such as e.g., displacement 
vectors Xi = \f\(t + At) — fi(t)\ between the 
time t and t + At for some fixed time interval 
At. We can also choose to be the vector 
displacement itself Xi = fi(t + At) — fi(t). 
If there are no correlations between the mo- 
tions of particles, then the variations of the 
XiS will satisfy 



(6) 



where a[x] denotes the mean square devia- 
tions of the random number x, a[x] = ((x— < 
x >) 2 ). However, some correlations between 
the particle motions will increase cr[J2i Xi] or 
anti-correlations will decrease it. Following 
Doliwa and Heuer, we can define the dynamic 
cooperativity as 



i\j-coop Aj 

x = EMX~_ 



(7) 



In the case of no correlations between the mo- 
tions of particles, as in (6), we get Nx° p = 1. 
If there exist some positive correlated mo- 
tions between particles, we would get A^ op > 
1, while anti-correlations between the mo- 
tions of particles would result in Nx° p < 1. 
Doliwa and Heuer investigated the dynamic 
cooperativity of hard sphere systems in 2D 
and 3D, where they found finite cooperativity 
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(Nx° p > 1) for both vector displacement and 
the scalar magnitude of the displacement, 
which is consistent with the snapshots of the 
particle motions in their work. They argue 
that the dynamic cooperativity measures the 
total reduction of degrees of freedom due to 
the correlations. Here we also studied the 
dynamic cooperativity of the lattice gas par- 
ticles by calculating iV£ >op for both the scalar 
displacement and the vector displacement it- 
self. Interestingly, we found that the scalar 
displacement exhibited finite dynamic coop- 
erativity (Fig. 9a), while the vector displace- 
ment itself showed weak anti-correlations be- 
tween particles, as shown in Fig. 9b. In the 
case of scalar displacement, the cooperativ- 
ity increases at first as the time interval At 
increases and reaches its maximum near the 
a-relaxation time scale r. Then it decrease 
back to values around unity (corresponding 
to no correlations) at large At. 

Contrasting features of cooperativity for 
our LCG system and that for the hard sphere 
systems may be interpreted as follows. In the 
case of hard sphere systems near the glass 
transition, the packing density is very high 
and the inter-particle interaction is a short 
ranged one. Therefore, the local motions of 
particles in hard sphere systems are naturally 
highly correlated in both its direction and 
magnitude due to the continuity constraint of 
particles resulting in a large scale flow with 
directional correlations. 

In contrast, in the case of the LCG, the 
density of particles is relatively low (/ ~ 
0.38) as compared with the case of hard 
sphere systems near the glass transition. In 
addition to that, charge motions in the LCG 
is driven by thermal effect. From the snap- 
shots of charge configurations, we see that 
there exist locally mobile regions as well as 
locally immobile regions. Locally immobile 
regions consist of charge configurations that 
are close to the low temperature striped pat- 
terns. Mobile regions, however, consist of 



charges that are agitated in random direc- 
tions due to the thermal effect. Thus we do 
not observe positive dynamic cooperativity in 
vector displacement, but only the scalar dis- 
placement exhibits appreciable positive coop- 
erativity due to the local regions with high 
mobilities. Hence, heterogeneity still exists 
in our lattice Coulomb gas in terms of local 
mobility distribution, but unlike the case of 
hard sphere systems, there is no appreciable 
average local flow. 

Also, we may look into the nature of the 
equilibrium dynamics of the system in wave- 
vector space. Figure 10 shows the structure 
factor S(q) = (\p q \ 2 ) at equilibrium where 
Pq = Ejexpfig • rj\/N where q = ^n, n = 
1, 2, • • • 2/L. We see that the structure factor 
of our LCG shows some similarity to those 
of dense liquids with first peak corresponding 
roughly to the inverse of the average distance 
between charges. Due to the lattice nature of 
the LCG, the wave vector has cutoff value at 
qmax = 7r as in the figure. 

The diffusive properties of the system can 
be probed by calculating the incoherent scat- 
tering function (ISF) Fs(q, t) which is defined 
as in our model of LCG 

F s (q,t) = QTexpig. fc(f) - f$(0)])/JV o , 

3=1 

(8) 

where fj(t) denotes the position of j-th parti- 
cle on the lattice. Due to the discrete lattice 
nature of our model sytem, we need to con- 
sider the wave-vectors within the first Bril- 
louin zone q = ^n, n = 0,1,2,- • • L — 1. 
Figure 11 shows the g-dependence of F s (q,t) 
at temperature T = 0.033. We find that the 
long-time behavior of Fs(q, t) also can be fit- 
ted to stretched exponential form. For low 
q, the late time (3 exponents were close to 
one (pure exponential relaxation) but as q 
increases the exponents decreased down to 
P « 0.73 for q = 18 x 2tt/36, and T = 0.033 
(Fig. 12). As can be seen from the definition 
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of Fs(q,t), for gaussian distribution for the 
displacement vector Ar i; we would get 



F G {q,t) = (expzg[Ar]> = exp[-|-((Ar) 2 )]. 



(9) 



Figure 13 shows that the gaussian approxi- 
mation is quite good for low q. That is, for 
long distance diffusion, the distribution gets 
closer to gaussian. However, as q becomes 
larger, the gaussian approximation gets worse 
as shown in the figure. Similar features were 
reported in molecular dynamics simulations 
on the dynamics of supercooled water |^5| . 

In summary, we have shown that the 2D 
LCG with fractional filling of charges exhibits 
an equilibrium relaxation behavior, above 
first order melting transition, characterized 
by two time-regimes of stretched exponetial 
form with temperature dependent exponents, 
which is quite similar to the primary re- 
laxation of typical supercooled liquids. We 
found a strong deviation from proportional- 
ity between the diffusive time scale and the 
structural relaxation time scale resembling 
the breakdown of SE relation in supercooled 
liquids. This is accompanied by a character- 
istic dynamic cooperativity, where the scalar 
displacement exhibits positive cooperativity 
while the vector displacement shows anti- 
correlations leading to the vector cooperativ- 
ity less than unity. We have identified the 
microscopic heterogeneous structure which is 
responsible for this phenomena. 
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FIGURE CAPTIONS 

Snapshots of charge configuration at 
time steps (a) t=16 MCS, (b) t=4096 
MCS, (c) t=65536 MCS, and (d) 
t=1048576 MCS, for temperature T = 
0.026 and / = 55/144, exhibiting 
coarsening toward an ordered striped 
state. Positive charges are represented 
by filled squares. 

Energy histogram near the first or- 
der transition temperature (for T = 
0.03165 and T = 0.0317). 

Regimes of charge patterns for the 
range of value of / between 1/3 and 
2/5. See the text for details. 

(a) The charge autocorrelation func- 
tions for temperatures T = 0.1, 0.08, 
0.06, 0.05, 0.042, 0.037, 0.035, 0.033. 

(b) Arrhenius plot for the relaxation 
time (log(r) versus 1/T). (c) Charge 
autocorrelation functions in (a) replot- 
ted in terms of the rescaled time t/r(T) 
which shows that the time-temperature 
superposition is broken, (d) Tempera- 
ture dependence of the b and f3 expo- 
nents for charge autocorrelation func- 
tions. 

Squared displacement W(t) versus time 
t for the same temperatures as in 
Fig. 4a. 

Comparison of the two time scales D^ 1 
and t (ADt versus T), which clearly 
shows that the diffusive time scale in- 
creases slowly (as the temperature is 
lowered) as compared with the struc- 
tural relaxation time. 

Typical charge configurations at T = 
0.033. Positive charges are represented 
by filled squares. 
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Fig. 8. Trajectories of moving positive charges 
at T = 0.033 over a time interval of 
500 MC steps. Arrows indicate the di- 
rections of single charge motions. 

Fig. 9. Dynamic cooperativity for (a) scalar 
displacement and (b) vector displace- 
ment respectively for varying time in- 
tervals at various temperatures. 

Fig. 10. The structure factor S{q) at T = 0.033 
and T = 0.037. 

Fig. 11. The incoherent intermediate scattering 
functions at temperature T = 0.033 for 
various wave vectors q. 

Fig. 12. g-dependence of the b and (3 exponents 
for the intermediate scattering func- 
tions at temperature T = 0.033. 

Fig. 13. Comparison of the Gaussian approxi- 
mations and the incoherent intermedi- 
ate scattering functions at temperature 
T = 0.033 for various wave vectors q. 
We can see that the gaussian approxi- 
mation is worse at large wave vectors. 
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(a) t=16MCS 




(b) t=4096MCS 




(c) t=65536MCS 




(d) t=1048576MCS 
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